# Author: Mark Richardson
# Purpose: Create Figure 2 - Predicted performance across the range of months vacant

# Load packages
library(dplyr)
library(estimatr)
library(ggplot2)
library(sjPlot)
library(patchwork)

# Load data
perf <- readxl::read_excel("data/01_lpr_012824.xlsx")

perf <- perf %>%
  mutate(mvac = (1+Total_Vacant_Days)/30,
         dpm = directPAS*mvac,
         emp = employ16/1000,
         lnemp = log(emp),
         mvsq = mvac^2)

# Model 2: Main model
model2 <- lm_robust(hier_perf_in ~ mvac + I(mvac^2) + directPAS + totturn + eop + cabinet +
                      offsec + indcom + depprior + skills_mean_2014 + ideo_rating,
                    data = perf, clusters= dept_acr, se_type = "stata")


summary(model2)


p1 <- plot_model(model2, type="pred", terms=c("mvac[all]"), colors = "black") +
  labs(x = "Months Vacant", y = "", title = "") +
  geom_rug(sides = "b", outside = TRUE) +
  theme_bw(base_size = 20) +
  scale_y_continuous(breaks = seq(from = -2, to = 1.5, by = 0.5),
                     limits = c(-2, 1.5)) +
  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(),
        axis.ticks.length.x = unit(0.4, "cm")) +
  coord_cartesian(clip = "off")

# Model 4: Direct PAS model
model4 <- lm_robust(hier_perf_in ~ mvac + I(mvac^2) + totturn + eop + cabinet +
                      offsec + indcom + depprior + skills_mean_2014 + ideo_rating,
                    data = perf %>% filter(directPAS == 1),
                    clusters= dept_acr, se_type = "stata")


summary(model4)

p2 <- plot_model(model4, type="pred", terms=c("mvac[all]"), colors = "black") +
  labs(x = "Months Vacant", y = "", title = "") +
  geom_rug(sides = "b", outside = TRUE) +
  theme_bw(base_size = 20) +
  scale_y_continuous(breaks = seq(from = -2, to = 1.5, by = 0.5),
                     limits = c(-2, 1.5)) +
  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(),
        axis.ticks.length.x = unit(0.4, "cm")) +
  coord_cartesian(clip = "off")


p1 + p2

ggsave("figs/fg2.png",
       device = "png", dpi = "retina",
       units = "in", width = 7, height = 5)
